/*
 * EquilibriumProblemReduced.cpp
 *
 *  Created on: 27 Jul 2011
 *      Author: Allan
 */

#include "EquilibriumProblemReduced.h"

// GeoReact includes
#include "EquilibriumSolver.h"

EquilibriumProblemReduced::EquilibriumProblemReduced(const EquilibriumSolver& eqSolver, double T, double P, const VectorXd& n) :
equilibriumSolver(eqSolver), jIndexes(eqSolver.jIndexes), eIndexes(eqSolver.eIndexes), vej(eqSolver.vej), alpha(eqSolver.alpha), beta(eqSolver.beta), 
Nj(eqSolver.Nj), Ne(eqSolver.Ne), T(T), P(P), n(n), a(n)
{
	// Allocate all the auxiliary Eigen objects with zero entries
	Ke = VectorXd::Zero(Ne);	Qe = VectorXd::Zero(Ne);
	ne = VectorXd::Zero(Ne);	h  = VectorXd::Zero(Nj);
	
	// Initialise the equilibrium constants of the equilibrium reactions (Ke)
	Ke = eqSolver.equilibriumReactions.EquilibriumConstants(T, P);
	
	// Initialise the vector (ne) with the values from the initial guess (n)
	for(uint e = 0; e < Ne; ++e) ne[e] = n[eIndexes[e]];
}

void EquilibriumProblemReduced::Function(const VectorXd& nj, VectorXd& F)
{
	// Update the auxiliary data members of the equilibrium problem
	UpdateAuxiliaryData(nj);
	
	// Compute the residual vector (F)
	F = alpha * nj + beta * ne - h;
}

void EquilibriumProblemReduced::Jacobian(const VectorXd& nj, MatrixXd& J)
{
	// Compute the jacobian matrix (J)
	for(uint i = 0; i < Nj; ++i) for(uint j = 0; j < Nj; ++j)
	{
		double sum = 0.0; for(uint e = 0; e < Ne; ++e) sum += beta(i, e) * vej(e, j) * ne[e];
		
		J(i, j) = alpha(i, j) + sum/nj[j];
	}
}

void EquilibriumProblemReduced::UpdateAuxiliaryData(const VectorXd& nj)
{
	// Check if all the entries in (nj) are non-negative
	if(nj.minCoeff() > 0)
	{
		// Update the nj part of the vector (n) with the values of the vector (nj)
		for(uint j = 0; j < Nj; ++j) n[jIndexes[j]] = nj[j];
		
		// Update the activity vector (a) with temperature (T), pressure (P) and the mole vector (n)
		a = equilibriumSolver.multiphase.Activities(T, P, n);
	
		// Update the reaction quotients of the equilibrium reactions (Qe)
		Qe = equilibriumSolver.equilibriumReactions.ReactionQuotients(a);
	
		// Update the moles of the equilibrium species (ne)
		ne = ne.array() * Qe.array() / Ke.array();
	
		// Update the ne part of the vector (n)
		for(uint e = 0; e < Ne; ++e) n[eIndexes[e]] = ne[e];
		
		// Update the vector (h), the right-hand side values of the speciation conditions
		for(uint j = 0; j < Nj; ++j) h[j] = equilibriumSolver.hfunctions[j](T, P, n, a);
	}
}
